library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.2     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(broom)
library(broom.mixed)
library(stats)
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
PUF_ELIX_IP2009 <- read.csv("PUF_ELIX_IP2009.csv")

PUF_ELIX_IP2009 <- within(PUF_ELIX_IP2009, {
  STATE <- factor(STATE)
  COUNTY <- factor(COUNTY)
  FULL_FIPS_CODE <- factor(FULL_FIPS_CODE)
  CBSA <- factor(CBSA)
})

PUF_ELIX_IP2009 <- PUF_ELIX_IP2009 %>%
  mutate(VISITS = TOTAL_VISITS>0)
Beginning Covariates:
BENE_AGE The beneficiary’s approximate age, calculated from the last day of the year 2009
B_SEX The beneficiary’s sex - 0 for males, 1 for females
B_DIED Whether or not this beneficiary died in 2009 (0 for No, 1 for Yes)
B_DIVERSE If the person is of color (1) or white (0)
HF_PROP_RANK The proportional rank of the person’s county against all the counties in their state on Health Factors (RWJF, 2010)
HO_PROP_RANK The proportional rank of the person’s county against all the counties in their state on Health Outcomes (RWJF, 2010)
TOTCHRONIC The max number of chronic conditions the beneficiary had in 2009
MEAN_ELIX_SCORE The mean Elixhauser score based on the person’s admissions in 2009 (Might be artificially 0 if a person with chronic conditions had no admissions in 2009.)

m1 <- glm(VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data=PUF_ELIX_IP2009, family=binomial(link="logit"))
summary(m1)
## 
## Call:
## glm(formula = VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + 
##     HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, 
##     family = binomial(link = "logit"), data = PUF_ELIX_IP2009)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6333  -0.4555  -0.2246  -0.1804   3.0203  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -3.699990   0.067367 -54.923  < 2e-16 ***
## BENE_AGE        -0.005143   0.000798  -6.445 1.16e-10 ***
## B_SEX            0.005323   0.020665   0.258  0.79674    
## B_DIVERSE       -0.051365   0.028395  -1.809  0.07046 .  
## B_DIED           0.251596   0.078843   3.191  0.00142 ** 
## HF_PROP_RANK     0.057217   0.054299   1.054  0.29200    
## HO_PROP_RANK    -0.066118   0.057123  -1.157  0.24709    
## TOTCHRONIC       0.413353   0.004435  93.200  < 2e-16 ***
## MEAN_ELIX_SCORE  0.475669   0.005182  91.794  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 93906  on 110329  degrees of freedom
## Residual deviance: 64097  on 110321  degrees of freedom
##   (1852 observations deleted due to missingness)
## AIC: 64115
## 
## Number of Fisher Scoring iterations: 6

In this model,
Age, Total # of Chronic Conditions and Mean Elixhauser Score are significant at the 99% CI against the outcome of having any inpatient “VISITS”.
As age goes up, the likelihood of any visit goes down. At some point, people die.
As the # of chronic conditions & mean Elixhauser scores increase the likelihood of any inpatient visits increase.

Whether or not the person Died in 2009 is significant at the 95% CI.
People often die in the hospital, so dying and inpatient visits are significantly related.

Not being white (B_DIVERSE) is significant at the 90% CI.
The estimate is negative, so there is a negative relationship between being not white and having visits.
Inverse relationships in correlation matrices show that there is a negative relationship between being a person of color and living longer, also system utilization / going to the doctor regularly.

Now that we’ve looked at the outcome of any “VISITS”, we can look at the outcome of TOTAL_VISITS>0

PUF_ELIX_IP2009_GT0 <- subset(PUF_ELIX_IP2009, TOTAL_VISITS>0)

The original data set had 113,305 observations, and this subset has 16,908.
First, let’s try OLS against the outcome of TOTAL_VISITS > 0. Then,
We will try a Poisson model with log link for the observations.

# Could be worse ... but standardized residual plots look pretty wonky!!
m0 <- lm(TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data=PUF_ELIX_IP2009_GT0)
summary(m0)
## 
## Call:
## lm(formula = TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + 
##     HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, 
##     data = PUF_ELIX_IP2009_GT0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0379 -0.3597 -0.1659  0.2019  2.8953 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.7239201  0.0287808  25.153   <2e-16 ***
## BENE_AGE        -0.0002269  0.0003318  -0.684    0.494    
## B_SEX            0.0069197  0.0085893   0.806    0.420    
## B_DIVERSE        0.0189714  0.0118204   1.605    0.109    
## B_DIED           0.3027413  0.0310598   9.747   <2e-16 ***
## HF_PROP_RANK     0.0006657  0.0224746   0.030    0.976    
## HO_PROP_RANK    -0.0140451  0.0235501  -0.596    0.551    
## TOTCHRONIC       0.0870507  0.0020450  42.568   <2e-16 ***
## MEAN_ELIX_SCORE  0.0468164  0.0020604  22.722   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5467 on 16723 degrees of freedom
##   (175 observations deleted due to missingness)
## Multiple R-squared:  0.1414, Adjusted R-squared:  0.141 
## F-statistic: 344.2 on 8 and 16723 DF,  p-value: < 2.2e-16
plot(m0)

The OLS / Gaussian model is a poor fit.

m2 <- glm(TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data=PUF_ELIX_IP2009_GT0, family=poisson(link = "log"))
summary(m2)
## 
## Call:
## glm(formula = TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + 
##     HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, 
##     family = poisson(link = "log"), data = PUF_ELIX_IP2009_GT0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9085  -0.3100  -0.1467   0.1304   2.1249  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -0.2020799  0.0466156  -4.335 1.46e-05 ***
## BENE_AGE        -0.0001485  0.0005293  -0.281    0.779    
## B_SEX            0.0049213  0.0137592   0.358    0.721    
## B_DIVERSE        0.0141917  0.0188282   0.754    0.451    
## B_DIED           0.2043975  0.0447024   4.572 4.82e-06 ***
## HF_PROP_RANK     0.0002444  0.0359756   0.007    0.995    
## HO_PROP_RANK    -0.0099202  0.0376879  -0.263    0.792    
## TOTCHRONIC       0.0681299  0.0033237  20.498  < 2e-16 ***
## MEAN_ELIX_SCORE  0.0357517  0.0032438  11.021  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3669.4  on 16731  degrees of freedom
## Residual deviance: 3034.6  on 16723  degrees of freedom
##   (175 observations deleted due to missingness)
## AIC: 39373
## 
## Number of Fisher Scoring iterations: 4

Age is no longer significant in this model.
Death, Chronic Conditions and Elixhauser Score remain significant.

Let’s compare this Poisson model to a Negative Binomial model.

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
m3 <- glm.nb(TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data=PUF_ELIX_IP2009_GT0)
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
summary(m3)
## 
## Call:
## glm.nb(formula = TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + 
##     B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, 
##     data = PUF_ELIX_IP2009_GT0, init.theta = 118176.1914, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9085  -0.3100  -0.1467   0.1304   2.1249  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -0.2020797  0.0466160  -4.335 1.46e-05 ***
## BENE_AGE        -0.0001485  0.0005293  -0.281    0.779    
## B_SEX            0.0049212  0.0137593   0.358    0.721    
## B_DIVERSE        0.0141917  0.0188283   0.754    0.451    
## B_DIED           0.2043973  0.0447029   4.572 4.82e-06 ***
## HF_PROP_RANK     0.0002444  0.0359759   0.007    0.995    
## HO_PROP_RANK    -0.0099201  0.0376882  -0.263    0.792    
## TOTCHRONIC       0.0681299  0.0033237  20.498  < 2e-16 ***
## MEAN_ELIX_SCORE  0.0357517  0.0032438  11.021  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(118176.2) family taken to be 1)
## 
##     Null deviance: 3669.3  on 16731  degrees of freedom
## Residual deviance: 3034.5  on 16723  degrees of freedom
##   (175 observations deleted due to missingness)
## AIC: 39375
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  118176 
##           Std. Err.:  170345 
## Warning while fitting theta: iteration limit reached 
## 
##  2 x log-likelihood:  -39354.98

And finally, let’s fit the original data to a Zero-Inflated Negative Binomial model

library(pscl)
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
m4 <- zeroinfl(TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, 
               data = PUF_ELIX_IP2009, dist = "negbin")
summary(m4)
## 
## Call:
## zeroinfl(formula = TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + 
##     HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data = PUF_ELIX_IP2009, 
##     dist = "negbin")
## 
## Pearson residuals:
##      Min       1Q   Median       3Q      Max 
## -1.68015 -0.34172 -0.11074 -0.08544 30.89873 
## 
## Count model coefficients (negbin with log link):
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -2.2930409  0.0527404 -43.478  < 2e-16 ***
## BENE_AGE        -0.0018493  0.0005757  -3.212  0.00132 ** 
## B_SEX            0.0246383  0.0150825   1.634  0.10235    
## B_DIVERSE        0.0212150  0.0206331   1.028  0.30385    
## B_DIED           0.3635194  0.0489567   7.425 1.13e-13 ***
## HF_PROP_RANK     0.0378638  0.0394239   0.960  0.33684    
## HO_PROP_RANK    -0.0648401  0.0414014  -1.566  0.11732    
## TOTCHRONIC       0.2441284  0.0035428  68.908  < 2e-16 ***
## MEAN_ELIX_SCORE  0.1583228  0.0041000  38.615  < 2e-16 ***
## Log(theta)      14.0213433 18.6855791   0.750  0.45302    
## 
## Zero-inflation model coefficients (binomial with logit link):
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.591730   0.212105   7.504 6.17e-14 ***
## BENE_AGE         0.008951   0.002618   3.420 0.000627 ***
## B_SEX            0.102321   0.064578   1.584 0.113090    
## B_DIVERSE        0.211687   0.087673   2.415 0.015756 *  
## B_DIED           0.316981   0.241339   1.313 0.189038    
## HF_PROP_RANK     0.047620   0.169204   0.281 0.778377    
## HO_PROP_RANK    -0.153982   0.178558  -0.862 0.388485    
## TOTCHRONIC      -0.223664   0.015084 -14.828  < 2e-16 ***
## MEAN_ELIX_SCORE -2.191743   0.080080 -27.369  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Theta = 1228547.7199 
## Number of iterations in BFGS optimization: 34 
## Log-likelihood: -4.173e+04 on 19 Df

As you can see, we could have jumped straight to the Zero-Inflated model, as it contains the hurdle model!

Let’s look at some of the plots.

# Diagnostic and Residuals plot for m2: Poisson regression (over dispersion is evident)
plot(m2)

# Diagnostic and Residuals plot for m3: Negative Binomial regression (still over dispersed)
plot(m3)

# Negative Binomial
ggplot(data = PUF_ELIX_IP2009_GT0, aes(x = TOTCHRONIC, y = TOTAL_VISITS)) +
  geom_jitter(height=0.05, width = 0.1) + 
  stat_smooth(method = "glm.nb")
## `geom_smooth()` using formula 'y ~ x'
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

# Poisson
ggplot(data = PUF_ELIX_IP2009_GT0, aes(x = TOTCHRONIC, y = TOTAL_VISITS)) +
  geom_jitter(height=0.05, width = 0.1) + 
  stat_smooth(method = "glm", method.args = c(family="poisson"))
## `geom_smooth()` using formula 'y ~ x'

They’re basically identical.

# Negative Binomial
ggplot(data = PUF_ELIX_IP2009_GT0, aes(x = MEAN_ELIX_SCORE, y = TOTAL_VISITS)) +
  geom_jitter(height=0.05, width = 0.1) + 
  stat_smooth(method = "glm.nb")
## `geom_smooth()` using formula 'y ~ x'
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

# Poisson
ggplot(data = PUF_ELIX_IP2009_GT0, aes(x = MEAN_ELIX_SCORE, y = TOTAL_VISITS)) +
  geom_jitter(height=0.05, width = 0.1) + 
  stat_smooth(method = "glm", method.args = c(family="poisson"))
## `geom_smooth()` using formula 'y ~ x'

These are very similar.
Based on the dispersion handling, the zero-inflated model was the most accurate thus far,
but let’s look at some different covariates, using the Poisson model.
What else could be related to the total # of inpatient visits?

# Can't use TOTAL_PHYS or TOTAL_CODES due to their near-perfect correlation with the outcome variable (see python corr matrix)
m2b <- glm(TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE + CVRG_MOS + RX_CVRG_MOS + TOTAL_LOS, data=PUF_ELIX_IP2009_GT0, family=poisson(link = "log"))
summary(m2b)
## 
## Call:
## glm(formula = TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED + 
##     HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE + 
##     CVRG_MOS + RX_CVRG_MOS + TOTAL_LOS, family = poisson(link = "log"), 
##     data = PUF_ELIX_IP2009_GT0)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.47013  -0.22946  -0.09558   0.07130   1.95714  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -2.654e-01  9.852e-02  -2.694  0.00706 ** 
## BENE_AGE        -3.128e-05  5.262e-04  -0.059  0.95260    
## B_DIVERSE        4.418e-03  1.884e-02   0.235  0.81455    
## B_DIED           1.462e-01  4.476e-02   3.267  0.00109 ** 
## HF_PROP_RANK    -1.167e-03  3.600e-02  -0.032  0.97414    
## HO_PROP_RANK    -3.371e-03  3.774e-02  -0.089  0.92883    
## TOTCHRONIC       4.371e-02  3.450e-03  12.670  < 2e-16 ***
## MEAN_ELIX_SCORE  1.959e-02  3.317e-03   5.906 3.51e-09 ***
## CVRG_MOS         1.675e-03  7.327e-03   0.229  0.81915    
## RX_CVRG_MOS      1.290e-03  1.720e-03   0.750  0.45322    
## TOTAL_LOS        2.829e-02  8.923e-04  31.703  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3669.4  on 16731  degrees of freedom
## Residual deviance: 2156.8  on 16721  degrees of freedom
##   (175 observations deleted due to missingness)
## AIC: 38499
## 
## Number of Fisher Scoring iterations: 4

Now, let’s see what impact State and County or CBSA would have.

m2c <- glm(TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE + TOTAL_LOS + STATE, data=PUF_ELIX_IP2009_GT0, family=poisson(link = "log"))
summary(m2c)
## 
## Call:
## glm(formula = TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED + 
##     HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE + 
##     TOTAL_LOS + STATE, family = poisson(link = "log"), data = PUF_ELIX_IP2009_GT0)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.45431  -0.22923  -0.09585   0.07105   1.95996  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -2.229e-01  6.572e-02  -3.391 0.000697 ***
## BENE_AGE                  -4.833e-05  5.282e-04  -0.091 0.927104    
## B_DIVERSE                  3.524e-03  1.929e-02   0.183 0.855078    
## B_DIED                     1.458e-01  4.489e-02   3.247 0.001165 ** 
## HF_PROP_RANK               3.099e-03  3.833e-02   0.081 0.935568    
## HO_PROP_RANK              -7.839e-03  4.055e-02  -0.193 0.846697    
## TOTCHRONIC                 4.386e-02  3.437e-03  12.759  < 2e-16 ***
## MEAN_ELIX_SCORE            1.957e-02  3.322e-03   5.892 3.83e-09 ***
## TOTAL_LOS                  2.832e-02  8.951e-04  31.642  < 2e-16 ***
## STATEAlaska                3.373e-02  1.783e-01   0.189 0.849975    
## STATEArizona              -3.518e-02  7.256e-02  -0.485 0.627764    
## STATEArkansas             -3.369e-03  7.795e-02  -0.043 0.965530    
## STATECalifornia           -1.741e-03  5.368e-02  -0.032 0.974129    
## STATEColorado             -1.584e-02  7.701e-02  -0.206 0.837059    
## STATEConnecticut          -7.701e-04  8.125e-02  -0.009 0.992437    
## STATEDelaware             -8.329e-03  1.266e-01  -0.066 0.947555    
## STATEDistrict of Columbia -4.154e-04  1.610e-01  -0.003 0.997941    
## STATEFlorida              -8.418e-03  5.432e-02  -0.155 0.876838    
## STATEGeorgia               1.598e-03  6.204e-02   0.026 0.979450    
## STATEHawaii               -2.357e-03  1.261e-01  -0.019 0.985087    
## STATEIdaho                 6.570e-03  1.124e-01   0.058 0.953389    
## STATEIllinois             -9.086e-03  5.785e-02  -0.157 0.875207    
## STATEIndiana              -1.101e-02  6.505e-02  -0.169 0.865647    
## STATEIowa                 -1.995e-02  7.905e-02  -0.252 0.800749    
## STATEKansas               -6.933e-02  8.411e-02  -0.824 0.409816    
## STATEKentucky             -9.691e-03  6.953e-02  -0.139 0.889152    
## STATELouisiana            -1.039e-02  7.232e-02  -0.144 0.885782    
## STATEMaine                -7.726e-03  9.895e-02  -0.078 0.937767    
## STATEMaryland              7.252e-03  6.953e-02   0.104 0.916924    
## STATEMassachusetts        -2.879e-02  6.451e-02  -0.446 0.655398    
## STATEMichigan             -6.699e-03  5.911e-02  -0.113 0.909776    
## STATEMinnesota             9.549e-03  7.240e-02   0.132 0.895067    
## STATEMississippi          -2.220e-02  7.367e-02  -0.301 0.763129    
## STATEMissouri             -1.074e-02  6.586e-02  -0.163 0.870468    
## STATEMontana              -3.387e-03  1.322e-01  -0.026 0.979564    
## STATENebraska             -7.422e-02  9.288e-02  -0.799 0.424217    
## STATENevada               -5.976e-02  9.925e-02  -0.602 0.547103    
## STATENew Hampshire        -8.496e-02  1.093e-01  -0.778 0.436822    
## STATENew Jersey            6.841e-05  6.161e-02   0.001 0.999114    
## STATENew Mexico           -1.578e-02  1.064e-01  -0.148 0.882135    
## STATENew York             -1.869e-02  5.575e-02  -0.335 0.737502    
## STATENorth Carolina       -1.092e-02  6.091e-02  -0.179 0.857768    
## STATENorth Dakota          8.357e-02  1.305e-01   0.640 0.521923    
## STATEOhio                 -4.039e-03  5.844e-02  -0.069 0.944893    
## STATEOklahoma              3.336e-02  7.448e-02   0.448 0.654260    
## STATEOregon               -2.961e-02  8.336e-02  -0.355 0.722473    
## STATEPennsylvania          2.958e-03  5.747e-02   0.051 0.958948    
## STATERhode Island         -4.776e-02  1.255e-01  -0.381 0.703450    
## STATESouth Carolina       -2.662e-03  6.984e-02  -0.038 0.969596    
## STATESouth Dakota         -4.225e-02  1.297e-01  -0.326 0.744577    
## STATETennessee            -4.676e-02  6.523e-02  -0.717 0.473507    
## STATETexas                 1.154e-02  5.500e-02   0.210 0.833866    
## STATEUtah                 -8.078e-03  1.036e-01  -0.078 0.937859    
## STATEVermont              -5.687e-02  1.552e-01  -0.366 0.714037    
## STATEVirginia             -3.215e-02  6.444e-02  -0.499 0.617847    
## STATEWashington           -4.231e-02  6.926e-02  -0.611 0.541304    
## STATEWest Virginia         7.322e-02  8.493e-02   0.862 0.388607    
## STATEWisconsin            -5.365e-04  6.801e-02  -0.008 0.993706    
## STATEWyoming              -1.627e-02  1.568e-01  -0.104 0.917394    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3669.4  on 16731  degrees of freedom
## Residual deviance: 2148.5  on 16673  degrees of freedom
##   (175 observations deleted due to missingness)
## AIC: 38587
## 
## Number of Fisher Scoring iterations: 4

No significant states. Let’s try CBSAs.

m2e <- glm(TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE + RX_CVRG_MOS + TOTAL_LOS + CBSA, data=PUF_ELIX_IP2009_GT0, family=Gamma(link = "log"))
summary(m2e)
## 
## Call:
## glm(formula = TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED + 
##     HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE + 
##     RX_CVRG_MOS + TOTAL_LOS + CBSA, family = Gamma(link = "log"), 
##     data = PUF_ELIX_IP2009_GT0)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.79865  -0.19922  -0.08349   0.06562   1.40956  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.200e-01  1.796e-02 -12.251  < 2e-16 ***
## BENE_AGE         1.427e-05  1.915e-04   0.075  0.94058    
## B_DIVERSE        4.897e-03  6.934e-03   0.706  0.48004    
## B_DIED           1.374e-01  1.793e-02   7.662 1.93e-14 ***
## HF_PROP_RANK     1.245e-02  1.840e-02   0.677  0.49865    
## HO_PROP_RANK    -1.783e-02  1.865e-02  -0.957  0.33882    
## TOTCHRONIC       3.665e-02  1.225e-03  29.930  < 2e-16 ***
## MEAN_ELIX_SCORE  1.704e-02  1.207e-03  14.116  < 2e-16 ***
## RX_CVRG_MOS      8.316e-04  6.053e-04   1.374  0.16948    
## TOTAL_LOS        3.161e-02  4.130e-04  76.528  < 2e-16 ***
## CBSA10180        2.063e-01  8.716e-02   2.367  0.01793 *  
## CBSA10420        1.350e-02  4.864e-02   0.277  0.78141    
## CBSA10500       -4.525e-02  9.031e-02  -0.501  0.61633    
## CBSA10580       -3.526e-03  4.807e-02  -0.073  0.94153    
## CBSA10740        5.629e-03  5.110e-02   0.110  0.91230    
## CBSA10780        2.706e-02  9.439e-02   0.287  0.77438    
## CBSA10900       -4.473e-02  4.071e-02  -1.099  0.27181    
## CBSA11020        6.080e-02  9.922e-02   0.613  0.54002    
## CBSA11100       -1.261e-01  8.681e-02  -1.452  0.14649    
## CBSA11180       -1.334e-01  1.399e-01  -0.954  0.34026    
## CBSA11260        3.009e-02  7.622e-02   0.395  0.69298    
## CBSA11300       -8.264e-02  8.688e-02  -0.951  0.34147    
## CBSA11340       -5.405e-02  8.690e-02  -0.622  0.53396    
## CBSA11460        9.511e-02  7.404e-02   1.285  0.19897    
## CBSA11500       -1.443e-02  9.900e-02  -0.146  0.88408    
## CBSA11540        1.259e-02  1.398e-01   0.090  0.92826    
## CBSA11700       -6.832e-02  5.497e-02  -1.243  0.21397    
## CBSA12020        8.107e-02  1.399e-01   0.580  0.56222    
## CBSA12060        5.941e-03  2.287e-02   0.260  0.79499    
## CBSA12100        1.355e-02  9.042e-02   0.150  0.88087    
## CBSA12220        4.768e-03  1.805e-01   0.026  0.97893    
## CBSA12260        2.515e-03  5.164e-02   0.049  0.96116    
## CBSA12420        1.120e-01  4.208e-02   2.661  0.00780 ** 
## CBSA12540       -5.175e-02  6.163e-02  -0.840  0.40110    
## CBSA12580        3.042e-04  2.642e-02   0.012  0.99082    
## CBSA12620       -1.091e-01  9.034e-02  -1.208  0.22720    
## CBSA12700        3.685e-02  6.683e-02   0.551  0.58131    
## CBSA12940        9.324e-03  5.406e-02   0.172  0.86307    
## CBSA12980        1.344e-01  1.106e-01   1.215  0.22440    
## CBSA13020       -2.486e-03  1.276e-01  -0.019  0.98446    
## CBSA13140       -2.607e-02  6.156e-02  -0.423  0.67196    
## CBSA13380       -1.036e-01  9.046e-02  -1.145  0.25228    
## CBSA13460        7.974e-02  1.044e-01   0.764  0.44496    
## CBSA13644        4.349e-02  4.572e-02   0.951  0.34146    
## CBSA13740       -1.727e-02  9.895e-02  -0.175  0.86143    
## CBSA13780        4.656e-02  7.193e-02   0.647  0.51748    
## CBSA13820        1.906e-02  3.535e-02   0.539  0.58968    
## CBSA13900        6.261e-03  9.901e-02   0.063  0.94958    
## CBSA13980       -6.874e-02  8.087e-02  -0.850  0.39533    
## CBSA14020       -1.059e-01  9.430e-02  -1.124  0.26123    
## CBSA14060        1.149e-02  1.183e-01   0.097  0.92260    
## CBSA14260       -4.319e-02  6.689e-02  -0.646  0.51847    
## CBSA14484       -2.021e-02  3.071e-02  -0.658  0.51044    
## CBSA14500        2.888e-02  9.050e-02   0.319  0.74961    
## CBSA14540        7.729e-02  1.398e-01   0.553  0.58039    
## CBSA14740       -1.168e-01  9.437e-02  -1.237  0.21602    
## CBSA14860       -4.225e-02  5.551e-02  -0.761  0.44665    
## CBSA15180        2.988e-02  6.651e-02   0.449  0.65327    
## CBSA15260        1.557e-01  1.563e-01   0.996  0.31928    
## CBSA15380        6.870e-02  4.076e-02   1.686  0.09191 .  
## CBSA15500       -8.876e-02  1.182e-01  -0.751  0.45263    
## CBSA15540       -2.147e-02  1.182e-01  -0.182  0.85587    
## CBSA15764       -2.191e-02  3.329e-02  -0.658  0.51046    
## CBSA15804        4.231e-04  3.828e-02   0.011  0.99118    
## CBSA15940        3.923e-03  5.925e-02   0.066  0.94721    
## CBSA15980        2.933e-03  4.748e-02   0.062  0.95074    
## CBSA16020        7.334e-03  8.369e-02   0.088  0.93017    
## CBSA16180        2.742e-01  1.562e-01   1.755  0.07921 .  
## CBSA16220        5.879e-02  1.105e-01   0.532  0.59483    
## CBSA16300       -6.216e-02  6.537e-02  -0.951  0.34169    
## CBSA16580       -8.000e-02  1.276e-01  -0.627  0.53079    
## CBSA16620        6.249e-04  6.561e-02   0.010  0.99240    
## CBSA16700        1.656e-01  5.183e-02   3.194  0.00140 ** 
## CBSA16740       -3.385e-02  3.427e-02  -0.988  0.32329    
## CBSA16820       -8.233e-02  7.835e-02  -1.051  0.29337    
## CBSA16860       -6.139e-02  5.332e-02  -1.151  0.24958    
## CBSA16940        3.225e-02  1.398e-01   0.231  0.81748    
## CBSA16974       -1.001e-02  1.621e-02  -0.617  0.53703    
## CBSA17020       -9.520e-03  1.182e-01  -0.081  0.93579    
## CBSA17140       -1.956e-02  2.827e-02  -0.692  0.48895    
## CBSA17300        1.421e-02  8.379e-02   0.170  0.86532    
## CBSA17420        1.277e-01  8.094e-02   1.578  0.11466    
## CBSA17460       -2.554e-02  2.778e-02  -0.920  0.35780    
## CBSA17660        1.720e-01  1.276e-01   1.348  0.17776    
## CBSA17780        2.605e-01  2.209e-01   1.179  0.23836    
## CBSA17820        3.677e-03  5.823e-02   0.063  0.94965    
## CBSA17860        3.178e-01  1.107e-01   2.871  0.00409 ** 
## CBSA17900       -2.428e-02  5.329e-02  -0.456  0.64862    
## CBSA17980       -1.088e-01  8.681e-02  -1.253  0.21014    
## CBSA18020        1.016e-01  1.106e-01   0.919  0.35834    
## CBSA18140        2.945e-02  3.284e-02   0.897  0.36986    
## CBSA18580        6.387e-02  5.924e-02   1.078  0.28099    
## CBSA18700        3.088e-02  3.124e-01   0.099  0.92126    
## CBSA18880        3.411e-02  1.044e-01   0.327  0.74386    
## CBSA19060       -3.655e-02  1.277e-01  -0.286  0.77464    
## CBSA19124        3.672e-02  2.629e-02   1.397  0.16254    
## CBSA19140       -1.809e-01  1.107e-01  -1.634  0.10224    
## CBSA19180       -2.675e-01  1.277e-01  -2.096  0.03614 *  
## CBSA19260        1.383e-01  9.040e-02   1.530  0.12592    
## CBSA19340        1.692e-01  7.008e-02   2.415  0.01575 *  
## CBSA19380       -1.943e-03  4.323e-02  -0.045  0.96415    
## CBSA19460       -8.210e-02  8.377e-02  -0.980  0.32710    
## CBSA19500       -6.733e-02  1.277e-01  -0.527  0.59792    
## CBSA19660       -2.127e-02  4.791e-02  -0.444  0.65703    
## CBSA19740       -5.435e-03  3.116e-02  -0.174  0.86153    
## CBSA19780       -1.536e-02  5.567e-02  -0.276  0.78266    
## CBSA19804       -4.410e-02  2.801e-02  -1.574  0.11545    
## CBSA20020        6.476e-02  8.688e-02   0.745  0.45608    
## CBSA20100       -1.171e-01  9.912e-02  -1.181  0.23763    
## CBSA20220        5.485e-02  1.182e-01   0.464  0.64273    
## CBSA20260       -1.287e-01  8.375e-02  -1.537  0.12426    
## CBSA20500       -1.839e-02  6.170e-02  -0.298  0.76567    
## CBSA20740       -1.299e-01  1.182e-01  -1.099  0.27167    
## CBSA20764        1.729e-02  2.808e-02   0.616  0.53808    
## CBSA20940        2.337e-01  1.107e-01   2.112  0.03472 *  
## CBSA21060        1.673e-01  1.107e-01   1.512  0.13061    
## CBSA21140        3.301e-02  1.048e-01   0.315  0.75272    
## CBSA21300       -1.290e-01  1.563e-01  -0.825  0.40928    
## CBSA21340       -6.132e-03  5.375e-02  -0.114  0.90916    
## CBSA21500        4.269e-02  7.637e-02   0.559  0.57614    
## CBSA21660       -1.452e-01  8.079e-02  -1.797  0.07241 .  
## CBSA21780       -5.174e-02  7.029e-02  -0.736  0.46162    
## CBSA21820       -2.118e-01  1.804e-01  -1.174  0.24044    
## CBSA22020       -1.231e-01  1.044e-01  -1.179  0.23830    
## CBSA22140       -6.267e-02  1.805e-01  -0.347  0.72849    
## CBSA22180        7.547e-02  7.184e-02   1.051  0.29348    
## CBSA22220        1.637e-02  6.703e-02   0.244  0.80706    
## CBSA22380       -8.428e-02  1.106e-01  -0.762  0.44606    
## CBSA22420       -3.014e-03  6.415e-02  -0.047  0.96252    
## CBSA22500       -7.560e-03  8.684e-02  -0.087  0.93062    
## CBSA22520       -1.068e-01  9.896e-02  -1.079  0.28066    
## CBSA22540        4.562e-02  1.276e-01   0.357  0.72077    
## CBSA22660       -3.477e-02  6.414e-02  -0.542  0.58778    
## CBSA22744       -1.128e-02  3.268e-02  -0.345  0.72999    
## CBSA22900        4.958e-02  6.842e-02   0.725  0.46866    
## CBSA23060        8.054e-02  6.151e-02   1.309  0.19040    
## CBSA23104        4.179e-03  3.533e-02   0.118  0.90584    
## CBSA23420        4.717e-02  5.655e-02   0.834  0.40418    
## CBSA23460       -7.451e-02  1.042e-01  -0.715  0.47471    
## CBSA23540        1.074e-01  8.685e-02   1.237  0.21622    
## CBSA23580        2.630e-01  9.448e-02   2.783  0.00539 ** 
## CBSA23844        2.916e-02  4.746e-02   0.614  0.53893    
## CBSA24020        7.196e-02  8.364e-02   0.860  0.38960    
## CBSA24140        6.842e-02  1.276e-01   0.536  0.59185    
## CBSA24220       -1.592e-01  1.277e-01  -1.247  0.21248    
## CBSA24300        7.302e-02  9.032e-02   0.808  0.41883    
## CBSA24340       -2.053e-02  4.374e-02  -0.469  0.63882    
## CBSA24500        1.412e-01  1.106e-01   1.277  0.20168    
## CBSA24540       -8.998e-02  1.563e-01  -0.576  0.56479    
## CBSA24580       -6.194e-02  9.034e-02  -0.686  0.49296    
## CBSA24660       -7.857e-02  4.983e-02  -1.577  0.11486    
## CBSA24780        1.647e-01  1.563e-01   1.054  0.29194    
## CBSA24860       -1.914e-02  4.714e-02  -0.406  0.68475    
## CBSA25060       -1.813e-01  6.838e-02  -2.652  0.00802 ** 
## CBSA25180        1.432e-01  6.538e-02   2.190  0.02852 *  
## CBSA25260       -2.269e-01  3.124e-01  -0.726  0.46767    
## CBSA25420        9.453e-03  5.390e-02   0.175  0.86080    
## CBSA25500        9.008e-02  9.440e-02   0.954  0.33999    
## CBSA25540        3.946e-02  3.796e-02   1.040  0.29855    
## CBSA25620       -4.831e-02  8.685e-02  -0.556  0.57804    
## CBSA25860        5.716e-02  5.388e-02   1.061  0.28871    
## CBSA25980       -3.831e-01  3.125e-01  -1.226  0.22017    
## CBSA26100       -3.510e-02  6.706e-02  -0.523  0.60067    
## CBSA26180       -6.240e-04  4.663e-02  -0.013  0.98932    
## CBSA26300       -4.481e-02  8.392e-02  -0.534  0.59343    
## CBSA26380       -6.305e-03  8.091e-02  -0.078  0.93789    
## CBSA26420        2.744e-02  2.275e-02   1.206  0.22789    
## CBSA26580       -3.036e-02  5.949e-02  -0.510  0.60985    
## CBSA26620        4.109e-02  6.706e-02   0.613  0.54003    
## CBSA26820       -7.029e-02  1.398e-01  -0.503  0.61515    
## CBSA26900        2.607e-02  3.136e-02   0.831  0.40579    
## CBSA26980        3.622e-01  1.805e-01   2.007  0.04481 *  
## CBSA27060       -1.592e-01  1.563e-01  -1.018  0.30847    
## CBSA27100       -9.886e-02  1.106e-01  -0.894  0.37135    
## CBSA27140        1.808e-02  5.661e-02   0.319  0.74937    
## CBSA27180       -6.696e-02  1.043e-01  -0.642  0.52081    
## CBSA27260       -3.862e-02  3.944e-02  -0.979  0.32751    
## CBSA27340       -1.816e-01  1.278e-01  -1.421  0.15531    
## CBSA27500        9.930e-03  1.276e-01   0.078  0.93799    
## CBSA27620        1.876e-02  1.106e-01   0.170  0.86528    
## CBSA27740        4.750e-02  8.377e-02   0.567  0.57075    
## CBSA27780       -1.093e-01  7.394e-02  -1.478  0.13948    
## CBSA27860       -1.800e-02  2.209e-01  -0.081  0.93505    
## CBSA27900       -8.769e-02  1.182e-01  -0.742  0.45823    
## CBSA28020       -2.265e-02  6.548e-02  -0.346  0.72940    
## CBSA28100        6.461e-02  1.182e-01   0.547  0.58461    
## CBSA28140        1.565e-03  2.962e-02   0.053  0.95787    
## CBSA28420        3.091e-02  8.366e-02   0.370  0.71174    
## CBSA28660        7.286e-02  7.604e-02   0.958  0.33799    
## CBSA28700        2.148e-02  5.935e-02   0.362  0.71746    
## CBSA28740        2.066e-01  1.105e-01   1.869  0.06166 .  
## CBSA28940       -3.699e-02  4.993e-02  -0.741  0.45887    
## CBSA29020       -1.136e-01  1.105e-01  -1.028  0.30416    
## CBSA29100        4.661e-02  9.055e-02   0.515  0.60673    
## CBSA29140        7.764e-02  1.043e-01   0.744  0.45673    
## CBSA29180        5.036e-02  8.380e-02   0.601  0.54785    
## CBSA29340        3.252e-01  1.563e-01   2.081  0.03743 *  
## CBSA29404        6.370e-02  4.757e-02   1.339  0.18055    
## CBSA29420        8.849e-03  8.689e-02   0.102  0.91889    
## CBSA29460        3.109e-02  4.596e-02   0.676  0.49875    
## CBSA29540       -8.085e-02  5.748e-02  -1.407  0.15955    
## CBSA29620       -6.053e-02  6.273e-02  -0.965  0.33462    
## CBSA29700       -3.211e-02  9.150e-02  -0.351  0.72565    
## CBSA29740        7.462e-03  9.921e-02   0.075  0.94004    
## CBSA29820       -3.746e-02  3.888e-02  -0.964  0.33522    
## CBSA29940       -2.719e-01  1.399e-01  -1.943  0.05202 .  
## CBSA30020        9.360e-02  1.107e-01   0.846  0.39774    
## CBSA30140        2.702e-02  1.277e-01   0.212  0.83236    
## CBSA30300       -2.240e-01  1.804e-01  -1.242  0.21438    
## CBSA30340        2.327e-01  1.398e-01   1.665  0.09587 .  
## CBSA30460       -3.149e-02  5.949e-02  -0.529  0.59664    
## CBSA30620       -1.250e-01  1.276e-01  -0.979  0.32738    
## CBSA30700       -2.356e-02  7.026e-02  -0.335  0.73739    
## CBSA30780        4.368e-02  5.480e-02   0.797  0.42546    
## CBSA30860       -1.408e-01  2.209e-01  -0.637  0.52408    
## CBSA30980       -2.462e-02  7.831e-02  -0.314  0.75325    
## CBSA31020       -1.641e-01  1.398e-01  -1.174  0.24056    
## CBSA31084        1.542e-02  1.912e-02   0.807  0.41993    
## CBSA31140       -2.318e-02  3.833e-02  -0.605  0.54527    
## CBSA31180        9.661e-02  1.043e-01   0.926  0.35446    
## CBSA31340       -7.304e-02  7.829e-02  -0.933  0.35086    
## CBSA31420        2.118e-02  1.106e-01   0.192  0.84812    
## CBSA31460        1.320e-01  1.398e-01   0.944  0.34503    
## CBSA31540       -6.989e-02  5.846e-02  -1.196  0.23190    
## CBSA31700       -8.889e-02  7.005e-02  -1.269  0.20447    
## CBSA31740       -2.770e-01  1.804e-01  -1.536  0.12466    
## CBSA31860       -4.094e-02  1.277e-01  -0.321  0.74845    
## CBSA31900       -2.156e-01  1.105e-01  -1.951  0.05110 .  
## CBSA32580       -7.456e-02  5.577e-02  -1.337  0.18124    
## CBSA32780        2.421e-01  1.182e-01   2.049  0.04045 *  
## CBSA32820        1.774e-02  3.939e-02   0.450  0.65246    
## CBSA32900        4.059e-02  1.563e-01   0.260  0.79507    
## CBSA33124       -3.579e-02  2.948e-02  -1.214  0.22474    
## CBSA33140       -9.469e-02  1.398e-01  -0.677  0.49816    
## CBSA33260       -2.298e-01  1.398e-01  -1.643  0.10042    
## CBSA33340        1.244e-02  3.229e-02   0.385  0.70017    
## CBSA33460       -7.118e-03  2.790e-02  -0.255  0.79859    
## CBSA33540       -3.547e-02  1.183e-01  -0.300  0.76424    
## CBSA33660       -3.726e-02  6.554e-02  -0.569  0.56970    
## CBSA33700        9.812e-03  6.537e-02   0.150  0.88069    
## CBSA33740       -1.647e-01  9.433e-02  -1.746  0.08077 .  
## CBSA33780        3.138e-01  1.106e-01   2.838  0.00455 ** 
## CBSA33860       -6.062e-02  6.155e-02  -0.985  0.32470    
## CBSA34060       -9.938e-02  1.278e-01  -0.778  0.43673    
## CBSA34100        1.771e-02  9.042e-02   0.196  0.84470    
## CBSA34580        1.653e-01  9.891e-02   1.671  0.09473 .  
## CBSA34620       -8.670e-02  1.043e-01  -0.831  0.40576    
## CBSA34740       -2.334e-02  8.369e-02  -0.279  0.78029    
## CBSA34820       -7.887e-02  7.597e-02  -1.038  0.29921    
## CBSA34900       -1.272e-02  3.123e-01  -0.041  0.96752    
## CBSA34940        8.095e-02  5.415e-02   1.495  0.13496    
## CBSA34980       -1.845e-02  3.331e-02  -0.554  0.57963    
## CBSA35004       -3.200e-02  2.602e-02  -1.230  0.21886    
## CBSA35084       -6.104e-03  2.809e-02  -0.217  0.82800    
## CBSA35300        2.518e-03  4.750e-02   0.053  0.95773    
## CBSA35380       -3.278e-02  4.079e-02  -0.804  0.42160    
## CBSA35644       -1.822e-02  1.440e-02  -1.265  0.20591    
## CBSA35660        1.230e-01  9.431e-02   1.304  0.19227    
## CBSA35840       -4.103e-02  3.678e-02  -1.116  0.26452    
## CBSA35980       -1.940e-03  7.824e-02  -0.025  0.98022    
## CBSA36084       -3.908e-02  3.370e-02  -1.160  0.24624    
## CBSA36100       -2.020e-02  6.272e-02  -0.322  0.74734    
## CBSA36140       -7.288e-02  9.436e-02  -0.772  0.43990    
## CBSA36220       -4.370e-02  1.182e-01  -0.370  0.71161    
## CBSA36260        4.754e-02  7.384e-02   0.644  0.51973    
## CBSA36420        5.650e-02  3.809e-02   1.483  0.13800    
## CBSA36500        1.404e-01  9.440e-02   1.488  0.13690    
## CBSA36540       -4.502e-02  4.687e-02  -0.961  0.33680    
## CBSA36740       -2.664e-02  2.904e-02  -0.917  0.35901    
## CBSA36780       -1.353e-01  9.454e-02  -1.431  0.15231    
## CBSA36980       -8.655e-03  9.440e-02  -0.092  0.92695    
## CBSA37100       -1.050e-02  5.104e-02  -0.206  0.83703    
## CBSA37340        5.896e-02  4.262e-02   1.383  0.16659    
## CBSA37380        1.774e-01  1.182e-01   1.501  0.13348    
## CBSA37460       -1.550e-01  9.435e-02  -1.642  0.10052    
## CBSA37620        1.127e-02  9.436e-02   0.119  0.90489    
## CBSA37700        3.985e-02  9.437e-02   0.422  0.67285    
## CBSA37764       -4.762e-02  4.455e-02  -1.069  0.28519    
## CBSA37860       -1.434e-02  4.965e-02  -0.289  0.77279    
## CBSA37900        9.432e-02  5.567e-02   1.694  0.09023 .  
## CBSA37964        2.685e-02  2.278e-02   1.178  0.23869    
## CBSA38060       -2.969e-02  2.604e-02  -1.140  0.25431    
## CBSA38220        2.313e-01  1.563e-01   1.480  0.13888    
## CBSA38300        2.492e-02  2.521e-02   0.989  0.32287    
## CBSA38340        2.879e-03  9.898e-02   0.029  0.97680    
## CBSA38540       -7.889e-02  1.277e-01  -0.618  0.53673    
## CBSA38860        1.769e-02  5.329e-02   0.332  0.73995    
## CBSA38900       -3.140e-02  3.183e-02  -0.986  0.32396    
## CBSA38940       -6.914e-02  5.557e-02  -1.244  0.21340    
## CBSA39100        8.094e-02  5.475e-02   1.478  0.13930    
## CBSA39140        5.048e-03  8.702e-02   0.058  0.95374    
## CBSA39300        2.622e-03  3.165e-02   0.083  0.93398    
## CBSA39340       -1.796e-01  9.442e-02  -1.903  0.05712 .  
## CBSA39380       -1.146e-01  1.106e-01  -1.037  0.29998    
## CBSA39460       -1.190e-01  8.369e-02  -1.422  0.15517    
## CBSA39540       -9.359e-02  8.088e-02  -1.157  0.24723    
## CBSA39580        4.262e-02  5.482e-02   0.777  0.43694    
## CBSA39660       -6.808e-03  9.431e-02  -0.072  0.94246    
## CBSA39740       -2.828e-03  7.013e-02  -0.040  0.96784    
## CBSA39820        3.652e-02  1.105e-01   0.330  0.74108    
## CBSA39900       -9.573e-02  7.015e-02  -1.365  0.17242    
## CBSA40060       -6.924e-02  4.002e-02  -1.730  0.08364 .  
## CBSA40140        1.062e-02  2.831e-02   0.375  0.70757    
## CBSA40220       -4.389e-02  7.003e-02  -0.627  0.53080    
## CBSA40340       -1.509e-01  1.277e-01  -1.182  0.23722    
## CBSA40380        3.795e-02  3.755e-02   1.011  0.31223    
## CBSA40420       -5.278e-03  6.401e-02  -0.082  0.93429    
## CBSA40484       -5.887e-02  6.036e-02  -0.975  0.32943    
## CBSA40580       -9.393e-02  1.562e-01  -0.601  0.54772    
## CBSA40660       -3.809e-02  9.896e-02  -0.385  0.70032    
## CBSA40900        4.263e-02  3.393e-02   1.257  0.20889    
## CBSA40980       -6.363e-02  9.440e-02  -0.674  0.50026    
## CBSA41060       -3.540e-03  1.043e-01  -0.034  0.97292    
## CBSA41100        7.611e-03  1.398e-01   0.054  0.95658    
## CBSA41140        3.651e-02  1.043e-01   0.350  0.72635    
## CBSA41180       -2.804e-02  2.482e-02  -1.130  0.25859    
## CBSA41420        1.483e-01  9.042e-02   1.641  0.10092    
## CBSA41500        1.055e-01  8.690e-02   1.215  0.22454    
## CBSA41540       -1.738e-01  1.106e-01  -1.571  0.11610    
## CBSA41620        6.831e-03  4.291e-02   0.159  0.87351    
## CBSA41660       -1.954e-02  9.038e-02  -0.216  0.82883    
## CBSA41700       -4.552e-03  3.325e-02  -0.137  0.89112    
## CBSA41740       -4.266e-02  2.977e-02  -1.433  0.15186    
## CBSA41780        6.907e-02  1.805e-01   0.383  0.70199    
## CBSA41884       -4.202e-02  3.204e-02  -1.312  0.18966    
## CBSA41940        7.531e-02  4.247e-02   1.774  0.07616 .  
## CBSA42020       -8.674e-02  8.376e-02  -1.036  0.30041    
## CBSA42044        1.263e-02  2.939e-02   0.430  0.66737    
## CBSA42060       -2.051e-02  5.931e-02  -0.346  0.72948    
## CBSA42100       -6.000e-02  9.904e-02  -0.606  0.54465    
## CBSA42140       -1.877e-01  1.277e-01  -1.469  0.14175    
## CBSA42220        3.460e-02  6.414e-02   0.539  0.58962    
## CBSA42340        6.356e-02  7.197e-02   0.883  0.37716    
## CBSA42540       -9.764e-02  4.624e-02  -2.112  0.03471 *  
## CBSA42644       -3.523e-02  3.018e-02  -1.167  0.24311    
## CBSA42680       -3.958e-02  9.045e-02  -0.438  0.66170    
## CBSA43100       -1.072e-01  1.182e-01  -0.907  0.36460    
## CBSA43300       -1.849e-01  1.563e-01  -1.183  0.23699    
## CBSA43340        1.269e-02  6.154e-02   0.206  0.83669    
## CBSA43580        1.943e-01  1.398e-01   1.390  0.16466    
## CBSA43620       -6.896e-02  8.097e-02  -0.852  0.39439    
## CBSA43780        8.605e-03  7.824e-02   0.110  0.91243    
## CBSA43900       -6.890e-02  7.388e-02  -0.933  0.35100    
## CBSA44060        3.548e-02  6.402e-02   0.554  0.57947    
## CBSA44100       -4.968e-02  7.832e-02  -0.634  0.52589    
## CBSA44140       -1.215e-02  4.182e-02  -0.291  0.77134    
## CBSA44180        8.930e-03  6.169e-02   0.145  0.88491    
## CBSA44220        2.007e-01  1.562e-01   1.285  0.19890    
## CBSA44300        2.067e-02  1.044e-01   0.198  0.84306    
## CBSA44600        2.170e-02  9.035e-02   0.240  0.81018    
## CBSA44700       -5.007e-02  6.851e-02  -0.731  0.46495    
## CBSA44940       -6.766e-02  1.043e-01  -0.649  0.51656    
## CBSA45060       -1.889e-02  5.239e-02  -0.361  0.71841    
## CBSA45104       -3.663e-02  5.093e-02  -0.719  0.47207    
## CBSA45220       -5.072e-02  7.830e-02  -0.648  0.51716    
## CBSA45300        4.195e-03  2.457e-02   0.171  0.86442    
## CBSA45460       -1.536e-01  8.682e-02  -1.769  0.07694 .  
## CBSA45500        7.446e-02  1.042e-01   0.714  0.47495    
## CBSA45780        1.436e-02  5.029e-02   0.286  0.77521    
## CBSA45820        1.336e-01  1.398e-01   0.956  0.33929    
## CBSA45940        6.435e-02  7.828e-02   0.822  0.41108    
## CBSA46060       -7.105e-02  4.226e-02  -1.681  0.09270 .  
## CBSA46140        2.594e-03  4.692e-02   0.055  0.95590    
## CBSA46220       -7.421e-02  7.391e-02  -1.004  0.31539    
## CBSA46340       -3.413e-02  7.598e-02  -0.449  0.65334    
## CBSA46540       -1.121e-02  8.702e-02  -0.129  0.89754    
## CBSA46660       -5.514e-02  1.182e-01  -0.467  0.64082    
## CBSA46700       -9.762e-02  7.825e-02  -1.248  0.21220    
## CBSA47020       -2.068e-01  1.277e-01  -1.619  0.10552    
## CBSA47220       -1.161e-01  1.183e-01  -0.982  0.32611    
## CBSA47260        8.331e-03  3.370e-02   0.247  0.80477    
## CBSA47300       -3.046e-02  9.043e-02  -0.337  0.73621    
## CBSA47380        9.507e-02  9.043e-02   1.051  0.29313    
## CBSA47580       -1.627e-01  1.183e-01  -1.376  0.16889    
## CBSA47644        7.739e-03  2.665e-02   0.290  0.77150    
## CBSA47894       -2.308e-02  2.555e-02  -0.903  0.36630    
## CBSA47940        3.127e-02  1.043e-01   0.300  0.76431    
## CBSA48140       -1.786e-01  1.398e-01  -1.277  0.20154    
## CBSA48300        1.892e-02  1.804e-01   0.105  0.91650    
## CBSA48424        3.792e-02  3.422e-02   1.108  0.26774    
## CBSA48540       -1.080e-02  9.031e-02  -0.120  0.90479    
## CBSA48620       -1.040e-01  5.478e-02  -1.898  0.05769 .  
## CBSA48660        3.701e-02  7.878e-02   0.470  0.63849    
## CBSA48700       -2.247e-02  1.182e-01  -0.190  0.84921    
## CBSA48864        8.832e-03  5.233e-02   0.169  0.86597    
## CBSA48900        6.705e-02  5.739e-02   1.168  0.24265    
## CBSA49020       -1.154e-01  1.398e-01  -0.826  0.40909    
## CBSA49180       -7.099e-02  6.843e-02  -1.037  0.29957    
## CBSA49340       -7.242e-02  5.389e-02  -1.344  0.17900    
## CBSA49420       -8.664e-02  9.894e-02  -0.876  0.38122    
## CBSA49620        8.453e-02  7.389e-02   1.144  0.25266    
## CBSA49660       -3.370e-02  5.313e-02  -0.634  0.52594    
## CBSA49700       -1.931e-01  1.563e-01  -1.236  0.21658    
## CBSA49740        2.856e-02  1.284e-01   0.222  0.82395    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.09748577)
## 
##     Null deviance: 2446.5  on 16731  degrees of freedom
## Residual deviance: 1342.3  on 16338  degrees of freedom
##   (175 observations deleted due to missingness)
## AIC: 12691
## 
## Number of Fisher Scoring iterations: 5

CBSA’s significant for the outcome of TOTAL_VISITS>0 include:

99% CI 12420 Austin-Round Rock, TX POSITIVE
16700 Charleston-North Charleston, SC POSITIVE
17860 Columbia, MO POSITIVE
23580 Gainesville, GA POSITIVE
25060 Gulfport-Biloxi-Pascagoula, MS NEGATIVE
95% CI 10180 Abilene, TX POSITIVE
19180 Danville, IL NEGATIVE
19340 Davenport-Moline-Rock Island, IA-IL POSITIVE
20940 El Centro, CA POSITIVE
25180 Hagerstown-Martinsburg, MD-WV POSITIVE
26980 Iowa City, IA POSITIVE
29340 Lake Charles, LA POSITIVE
90% CI 15380 Buffalo-Cheektowaga-Niagara Falls, NY POSITIVE
16180 Carson City, NV POSITIVE
21660 Eugene, OR NEGATIVE
28740 Kingston, NY POSITIVE

Now, let’s plot the outcome of TOTAL_VISITS against the MEAN_ELIX_SCORE, creating a separate intercept for each CBSA.

ggplot(data=PUF_ELIX_IP2009_GT0, aes(x=MEAN_ELIX_SCORE, y=TOTAL_VISITS, 
                                     group=CBSA)) + 
  geom_line() + 
  geom_smooth(method = "glm", method.args=c(family="poisson"), se=FALSE)
## `geom_smooth()` using formula 'y ~ x'

#
# Poisson against the outcome of TOTAL_LOS (Total Length of Stay) for each CBSA

ggplot(data=PUF_ELIX_IP2009_GT0, aes(x=MEAN_ELIX_SCORE, y=TOTAL_LOS, 
                                     group=CBSA)) + 
  geom_line() + 
  geom_smooth(method = "glm", method.args=c(family="poisson"), se=FALSE)
## `geom_smooth()` using formula 'y ~ x'

Now, let’s try fitting some models with the random effect of CBSA (noting that there are several beneficiaries whose counties are too small to belong to a CBSA)

zero_cbsas <- PUF_ELIX_IP2009 %>% subset(CBSA=="0")
zero_fips <- PUF_ELIX_IP2009 %>% subset(FULL_FIPS_CODE=="0")

25,943 beneficiaries without CBSAs.
1,755 beneficiaries without FIPS.

PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
  mutate(MC_SCALED = scale(MEAN_CHRONIC))

PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
  mutate(TC_SCALED = scale(TOTAL_CODES))

PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
  mutate(LOS_SCALED = scale(TOTAL_LOS))

PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
  mutate(PHYS_SCALED = scale(TOTAL_PHYS))

PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
  mutate(ELIX_SCALED = scale(MEAN_ELIX_SCORE))

CBSA_GT0 <- subset(PUF_ELIX_IP2009_GT0, CBSA!="0")

Modeling 13,080 Beneficiaries with 1 or more inpatient visits in 2009 who live in a CBSA:

# Linear Random Effects Model ~ RE Slope: Elixhauser Score, RE Intercept: CBSA
lmm = lmer(TOTAL_VISITS ~ MC_SCALED + LOS_SCALED + ELIX_SCALED + B_DIED + (ELIX_SCALED|CBSA), data = CBSA_GT0) 
## boundary (singular) fit: see ?isSingular
summary(lmm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: TOTAL_VISITS ~ MC_SCALED + LOS_SCALED + ELIX_SCALED + B_DIED +  
##     (ELIX_SCALED | CBSA)
##    Data: CBSA_GT0
## 
## REML criterion at convergence: 16922
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6815 -0.5634 -0.1767  0.2737  5.8067 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr
##  CBSA     (Intercept) 0.0002498 0.01580      
##           ELIX_SCALED 0.0002179 0.01476  1.00
##  Residual             0.2124769 0.46095      
## Number of obs: 13080, groups:  CBSA, 384
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 1.297e+00  4.292e-03 1.824e+02 302.116  < 2e-16 ***
## MC_SCALED   1.104e-01  4.192e-03 1.307e+04  26.331  < 2e-16 ***
## LOS_SCALED  3.032e-01  4.222e-03 1.307e+04  71.807  < 2e-16 ***
## ELIX_SCALED 5.054e-02  4.330e-03 2.062e+02  11.672  < 2e-16 ***
## B_DIED      2.188e-01  2.921e-02 1.307e+04   7.491 7.27e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) MC_SCA LOS_SC ELIX_S
## MC_SCALED   -0.001                     
## LOS_SCALED   0.004 -0.232              
## ELIX_SCALED  0.095 -0.087 -0.156       
## B_DIED      -0.132  0.002 -0.045  0.001
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
plot(lmm)

# GLM Random Effects Model, Gamma ~ RE Intercept: CBSA
glmm = glmer(TOTAL_VISITS ~ MC_SCALED + LOS_SCALED + ELIX_SCALED + B_DIED + (1|CBSA), data = CBSA_GT0, family = Gamma(link = "log")) 
summary(glmm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: TOTAL_VISITS ~ MC_SCALED + LOS_SCALED + ELIX_SCALED + B_DIED +  
##     (1 | CBSA)
##    Data: CBSA_GT0
## 
##      AIC      BIC   logLik deviance df.resid 
##   9713.2   9765.6  -4849.6   9699.2    13073 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8196 -0.5961 -0.2683  0.1865  6.5000 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CBSA     (Intercept) 0.001422 0.03771 
##  Residual             0.098513 0.31387 
## Number of obs: 13080, groups:  CBSA, 384
## 
## Fixed effects:
##             Estimate Std. Error t value Pr(>|z|)    
## (Intercept) 0.227911   0.004398  51.825  < 2e-16 ***
## MC_SCALED   0.077905   0.002604  29.915  < 2e-16 ***
## LOS_SCALED  0.193077   0.002766  69.793  < 2e-16 ***
## ELIX_SCALED 0.037110   0.002605  14.243  < 2e-16 ***
## B_DIED      0.132857   0.018120   7.332 2.27e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) MC_SCA LOS_SC ELIX_S
## MC_SCALED    0.002                     
## LOS_SCALED   0.003 -0.263              
## ELIX_SCALED  0.004 -0.104 -0.178       
## B_DIED      -0.079 -0.003 -0.046 -0.002
plot(glmm)

anova(lmm, glmm)
## refitting model(s) with ML (instead of REML)

In this case, the Gamma mixed model has the best fit, with AIC 9713.2 and -log-likelihood -4849.6
Both models have the same number of parameters.

CBSA_GT0 <- CBSA_GT0 %>%
  mutate(SCALED_COSTS = scale(ALLCOSTS))

# Random Effects Model ~ RE Slope: TOTAL_PHYS (Number of Physicians), RE Intercept: CBSA
lmm2 = lmer(TOTAL_LOS ~ MEAN_CHRONIC + TOTAL_CODES + TOTAL_PHYS + MEAN_ELIX_SCORE + SCALED_COSTS + (TOTAL_PHYS|CBSA), data = CBSA_GT0) 
## boundary (singular) fit: see ?isSingular
summary(lmm2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## TOTAL_LOS ~ MEAN_CHRONIC + TOTAL_CODES + TOTAL_PHYS + MEAN_ELIX_SCORE +  
##     SCALED_COSTS + (TOTAL_PHYS | CBSA)
##    Data: CBSA_GT0
## 
## REML criterion at convergence: 77595.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4172 -0.5529 -0.1889  0.3470  9.0341 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  CBSA     (Intercept)  4.1445  2.0358        
##           TOTAL_PHYS   0.9028  0.9502   -1.00
##  Residual             21.4792  4.6346        
## Number of obs: 13080, groups:  CBSA, 384
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     9.136e-01  2.170e-01 6.463e+02   4.211 2.91e-05 ***
## MEAN_CHRONIC    1.172e-01  2.097e-02 1.287e+04   5.590 2.31e-08 ***
## TOTAL_CODES     3.303e-01  1.472e-02 1.305e+04  22.447  < 2e-16 ***
## TOTAL_PHYS      2.405e-01  9.904e-02 8.459e+02   2.428   0.0154 *  
## MEAN_ELIX_SCORE 9.563e-02  2.058e-02 1.287e+04   4.646 3.42e-06 ***
## SCALED_COSTS    1.717e+00  7.058e-02 1.284e+04  24.330  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) MEAN_C TOTAL_C TOTAL_P MEAN_E
## MEAN_CHRONI -0.310                              
## TOTAL_CODES -0.064 -0.160                       
## TOTAL_PHYS  -0.659 -0.034 -0.535                
## MEAN_ELIX_S -0.181 -0.041 -0.180   0.048        
## SCALED_COST  0.516  0.057 -0.328  -0.233   0.023
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
plot(lmm2)

CBSA plus the slope of number of physicians explains the majority of the variance.
The remaining variance is small, but is well-explained by the covariates.
The residuals are fairly equally distributed around the zero line.
Let’s scale variables, and compare it with a glmer.

# GLM Random Effects Model ~ RE Slope: PHYS_SCALED (Number of Physicians), RE Intercept: CBSA
glmm2 = glmer.nb(TOTAL_LOS ~ MC_SCALED + TC_SCALED + PHYS_SCALED + ELIX_SCALED + SCALED_COSTS + (PHYS_SCALED|CBSA), data = CBSA_GT0) 
## boundary (singular) fit: see ?isSingular
summary(glmm2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(3.4561)  ( log )
## Formula: TOTAL_LOS ~ MC_SCALED + TC_SCALED + PHYS_SCALED + ELIX_SCALED +  
##     SCALED_COSTS + (PHYS_SCALED | CBSA)
##    Data: CBSA_GT0
## 
##      AIC      BIC   logLik deviance df.resid 
##  69576.5  69651.3 -34778.3  69556.5    13070 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6393 -0.7309 -0.2586  0.4263  8.5202 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev.  Corr 
##  CBSA   (Intercept) 7.041e-04 2.654e-02      
##         PHYS_SCALED 1.698e-10 1.303e-05 -0.87
## Number of obs: 13080, groups:  CBSA, 384
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.743218   0.006475 269.202  < 2e-16 ***
## MC_SCALED    0.049163   0.006550   7.505 6.12e-14 ***
## TC_SCALED    0.332325   0.013244  25.093  < 2e-16 ***
## PHYS_SCALED  0.003009   0.012166   0.247    0.805    
## ELIX_SCALED  0.040299   0.006256   6.442 1.18e-10 ***
## SCALED_COSTS 0.121043   0.009575  12.641  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) MC_SCA TC_SCA PHYS_S ELIX_S
## MC_SCALED   -0.019                            
## TC_SCALED   -0.052 -0.156                     
## PHYS_SCALED  0.006 -0.054 -0.677              
## ELIX_SCALED -0.013 -0.025 -0.182  0.059       
## SCALED_COST -0.005  0.056 -0.355 -0.262  0.024
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
plot(glmm2)

ANOVA Comparison:

anova(lmm2, glmm2)
## refitting model(s) with ML (instead of REML)

The glmer NB model has the superior fit, but this model doesn’t fit nearly as well as the visits model.

# LA COUNTY 06 037
county_subset <- CBSA_GT0 %>%
    mutate(predictions = predict(glmm)) %>%
      filter(FULL_FIPS_CODE=="6037") 

Exploring LA County’s 320 Beneficiaries, using the TOTAL_VISITS Gamma model

ggplot(county_subset, aes(x = predictions, y = TOTAL_VISITS)) + 
  geom_point() + 
  xlab("Predicted TOTAL VISITS (glmer Gamma)") + 
  ylab("Total Visits") + 
  theme_classic()

ggplot(county_subset, aes(x =predictions, y = (TOTAL_VISITS-predictions))) + 
  geom_point() + 
  xlab("Predicted Total Visits (glmer Gamma)") + 
  ylab("Residuals") + 
  theme_classic()

The best-fitting model still under-predicts in general

ggplot(county_subset, aes(x = MC_SCALED, y = TOTAL_VISITS)) + 
  geom_point() + 
    stat_smooth(method = "lm", aes(color="linear"), se=FALSE) +
    stat_smooth(method = "glm", aes(color="gamma"), method.args=list(family = Gamma(link = "log")), se=FALSE) +
    stat_smooth(method = "glm.nb", aes(color="negative binomial"), se=FALSE) +
  xlab("Mean Chronic Conditions (Scaled)") + 
  ylab("Number of Visits") + 
  theme_classic()
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

ggplot(county_subset, aes(x = ELIX_SCALED, y = TOTAL_VISITS)) + 
  geom_point() + 
    stat_smooth(method = "lm", aes(color="linear"), se=FALSE) +
    stat_smooth(method = "glm", aes(color="gamma"), method.args=list(family = Gamma(link = "log")), se=FALSE) +
    stat_smooth(method = "glm.nb", aes(color="negative binomial"), se=FALSE) +
  xlab("Avg Elix Score (Scaled)") + 
  ylab("Total Visits") + 
  theme_classic()
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

ggplot(county_subset, aes(x = LOS_SCALED, y = TOTAL_VISITS)) + 
  geom_point() + 
    stat_smooth(method = "lm", aes(color="linear"), se=FALSE) +
    stat_smooth(method = "glm", aes(color="gamma"), method.args=list(family = Gamma(link = "log")), se=FALSE) +
    stat_smooth(method = "glm.nb", aes(color="negative binomial"), se=FALSE) +
  xlab("Length of Stay (Scaled)") + 
  ylab("Total Visits") + 
  theme_classic()
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached